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Abstract 

An hyperelastic biphasic model is presented. For slow-draining 
problems (permeability less than lxl0~ 2 mm^V 1 ), numerical 
instabilities in the form of non-physical oscillations in the pressure field 
are observed in 3D problems using tetrahedral Taylor-Hood finite elements. 
As an alternative to considerable mesh refinement, a Galerkin least-square 
stabilization framework is proposed. This technique drastically reduces 
the pressure discrepancies and prevents these oscillations from propagating 
towards the centre of the medium. The performance and robustness of this 
technique are demonstrated on a 3D numerical example. 

Keywords: stabilization, porous media, soft tissue, Galerkin least-square, hyperelastic, 
large strains 

1 Introduction 

The theory of porous media (TPM) is a powerful and yet simple tool to model 
multi-phase media, primarily comprising a solid and a fluid (which are referred 
to as "constituents"). Due to its modularity, this theory also represents a practical 
framework to include other effects such as mass exchange, chemical reactions 
and electrochemical phenomena. There is no restriction as to which material 
models can be used, e.g. anisotropy or viscosity and is therefore well suited for 
the analysis of hydrogels, polymeric foams and hydrated biological soft tissues 
(e.g. cartilage and intervertebral disc). 

The concepts behind the TPM find their roots in diffusion and soil mechanics 
problems formulated in the nineteenth century (for historical developments, see 
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[3] and [9]). The TPM is a homogenized macroscopic representation of the porous 
media. It is a continuum-based model which fully couples the fluid and the solid 
(see Fig. 1). It overcomes the difficulty of obtaining an accurate geometrical 
description of the microstructure by using the concept of volume fractions to 
"smear" the constituent properties over a control space to obtain properties of 
the overall mixture. This is described in the following section. 
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Figure 1: The TPM representation (right) of the microstructure (left); at each 
material point, fluid and solid coexist. Displacement and pressure fields are 
coupled. 

The authors have employed TPM in the development of an hyperelastic biphasic 
swelling model for modelling the intervertebral disc. This application illustrates 
slow-draining problems, where a porous medium with low permeability (in the 
range 1 x 10~ 2 — 1 x 10~ 3 mm 4 N _1 s _1 ) is rapidly loaded. Such low permeabilities 
are typical for soft tissues (e.g. cartilage, brain tissue), mortar or homogeneous 
clays ([2], [13], [17]). 

This model has been implemented in a finite element framework, employing 
Taylor-Hood (quadratic shape functions for the solid displacement, linear shape 
functions for the pressure) tetrahedral elements, aiming to fulfil the inf-sup 
condition (see [6] and [5]). However, numerical instabilities manifest in the 
form of non-physical oscillations in the pressure field, which is a shortcoming 
already observed in the past (see for example [16] for a recent review of biological 
applications). Vermeer and Verruijt [20] explain that these instabilities occur 
because loads applied to free-flow boundary conditions may lead to singularities 
in the derivatives of the pressure field. They also derive a lower bound critical 
time-step for one-dimensional problems, suggesting the requirement for large 
time-steps to overcome this issue, often incompatible with fast loading rates. 

Several stabilisation techniques have been proposed in the context of Biot's 
consolidation problems for small deformations (e.g. [12] using least-squares 
mixed finite element methods, and [1] by perturbation of the flow equation). The 
current work proposes to stabilize the pressure oscillations in the context of TPM 
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for finite deformation problems, using a Galerkin Least-Square (GLS) formulation 
based on [19]. In order to focus only on the stabilization aspect, only a biphasic 
mixture is considered, that is a porous medium composed of two constituents a: 
an isotropic, non viscous hyperelastic solid (a = S) and an ideal fluid (a = F). 



2 The biphasic model 

The biphasic model presented in this section is based on [8] and [10]. A few 
preliminary assumptions are made to keep the derivation as simple as possible 
in order to focus on the stabilization aspect. First, the quasi-static problem 
of small biological tissues is herein considered, thus neglecting external body 
forces. Second, the constituents are assumed immiscible and no density supply 
is allowed. Third, it is assumed that the whole space is occupied by either 
of the constituents. Finally, intrinsic incompressibility is assumed for both 
constituents. It is important to note that the last assumption is not equivalent 
to incompressibility of the whole mixture. 

The most fundamental concept of mixture theory (established as early as in [15]) 
asserts that at any time t and at each spatial point x of the continuum, particles 
of both constituents a coexist. This implies that any elementary volume dv 
is simultaneously occupied by both phases and is split into partial elementary 
volumes dv a . The volume fractions can then be defined as: 

n°(x,t) = — a = {S,F} (2.1) 
Assuming that there is no gas in the mixture, the saturation condition is: 

n s + n F = 1 (2.2) 

The useful relationship between the apparent density p a and true density p% of 
each constituent is defined as follows: 

^ dv dv a dv ^ T 

The other main principles of mixture theory are summarised in "Truesdell's 
metaphysical principles" [18] l . The first consequence of these principles is that 



'"Truesdell's metaphysical principles": 1) All properties of the mixture must be mathematical 
consequences of properties of the constituents; 2) So as to describe the motion of a constituent, 
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each constituent a of the mixture is described by its own state of motion x a ( see 
Eq. 2.4), relating at time t the position vector X Q of a particle in the reference 
configuration to its position in the current configuration x Q . 

x* = f(X a ,i) (2.4) 

This subsequently implies the existence of a velocity field v a and a deformation 
gradient F a for each constituent. 

It proves convenient to describe the fluid velocity as a velocity relative to the solid 
constituent, leading to the introduction of the seepage velocity w: 

w = n F (v 5 - v F ) (2.5) 

As mentioned earlier, each constituent's intrinsic incompressibility does not imply 
mixture incompressibility. Any TPM elementary volume dv can be thought of as 
delimited by the porous solid; therefore, when subjected to deviatoric strains, the 
volume change of an elementary volume dv can directly be expressed as a function 
of the seepage velocity. 

Furthermore, the mixture as a whole can be seen as the superposition of its 
constituents. Therefore both the balance of mass (Eq. 2.6) and the balance of 
linear momentum (Eq. 2.7) are first expressed for each constituent independently. 
The rate at which momentum is transmitted by constituent a to the other 
constituent is accounted for in p a . 

^ + div (p a w a ) =0 a = {S, F} (2.6) 
at 

div (T a ) + p a = a = {S, F} (2.7) 

It is worth noting here that the apparent stress T a has the unit of a force per unit 
mixture area, not unit area of constituent a. 

The governing equations of both constituents are then superimposed in order to 
obtain the governing equation of the mixture as a whole. Using (Eq. 2.3) and the 
fact that the true density is constant (assumption of intrinsic incompressibility), 
the mass balance equation for the mixture can be written as: 

we may in imagination isolate it from the rest of the mixture, provided we allow properly for the 
actions of the other constituents upon it; 3) The motion of the mixture is governed by the same 
equations as is a single body. 
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— ^— ^ + div (n s v s + ti v ) = (2.8) 

ot v ' 

Using the saturation condition (Eq. 2.2) and the seepage velocity (Eq. 2.5) allows 
us to rewrite the mass balance equation of the mixture: 

div (v 5 + w) = (2.9) 

It is interesting to realise that as a result of the "smearing" introduced by the 
concept of volume fractions, the seepage velocity w represents a macro-level 
average of the fluid velocity with respect to the solid phase and not an actual 
fluid velocity at pore-level ([7]). 

In a similar fashion, the linear momentum balance of the mixture is obtained by 
summing (Eq. 2.7) for both constituents: 

div (T F + T s ) + p F + p s = (2. 10) 

The introduction of the volume fractions has the consequence of introducing one 
extra variable for each phase of the mixture, which implies that the "closure 
problem" is not fulfilled ([8]) (i.e. a greater number of unknowns than equations 
available). In order to overcome this problem, the saturation condition (Eq. 2.2) 
is imposed as a kinematic constraint onto the entropy inequality, introducing a 
Lagrange multiplier p (see for example [14] or [9] and references therein for 
complete derivation). As a direct consequence, the concept of effective stress 
is introduced: 

T a = -n a pl + T E (2.11) 

The quantities T E are the constituent effective stresses, which must be determined 
constitutively. Summing over constituents, the mixture stress T is expressed in 
terms of the mixture effective stress T E '- 



T = T + T = —pi + T E (2.12a) 

■ E 



T S = T| + T| (2.12b) 



It is common (see for example discussion in [1 1]) to consider the fluid as ideal (i.e. 
neglecting viscosity for slow draining materials such as soft tissues or hydrogels) 
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and neglect the dissipative stress in comparison to the drag generated by 
fluid-solid interactions p F , which is accounted for with Darcy's law. Hence, 
introducing <r as a new notation for the solid effective stress and identifying p 
as the fluid pressure, (Eq. 2.12a) simplifies to: 



T = -pI + cr (2.13) 

Truesdell's third metaphysical principle implies that the sum over all constituents 
of the partial balance relations has to take the same form as the balance of 
the single-phase material (i.e. in the current application divT = 0). As a 
consequence, the constituent momentum exchanges cancel out: p F + p s = and 
the linear momentum equation of the mixture is finally obtained by substituting 
(Eq. 2.13) into(Eq. 2.10): 



div {a- pi) =0 (2.14) 

The weak formulation of the problem is derived in a standard fashion (see Eq. 2. 15 
for expression in the current configuration), f and g are the weighting functions of 
the linear momentum and mass balance respectively. In (Eq. 2.15a) n represents 
the outward normal vector and in (Eq. 2.15b) t, is the surface traction vector. 



/ (g div (v) — w V<?) dv — gwnda (2.15a) 

J v J a 

f (Vf) : (tr-pl) dv= I f.tda (2.15b) 

J v J a 



Finally, the standard Galerkin formulation is obtained by discretising the weak 
form in space using Taylor-Hood elements, where a linear approximations for the 
pressure field (four-node tetrahedron) and a quadratic one for the displacement 
field (ten-node tetrahedron) are used. The weighting functions are discretised 
using the same shape functions as those used for displacement and pressure. The 
discretised quantities can be expressed in matrix form as: 

p = N p p e u = N"u e 

# = N p g e f = N u f e (2.16) 
The variational form of (Eq. 2.15) is: 

R„= f {NY} T VNVrfH I k {VN p g e } T {VN p p 6 } dv- 

J V J V 
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f {N p g e } T q^a = (2.17a) 

J a 

R u = I {VN u f e } T (<r e + N p p e ) dv - I {N"f e } T tck = (2.17b) 

J v J a 



R = R u + AtR p = (2.17c) 



A backward finite difference scheme is used for time integration (i.e. u t +At = 
(ut+At — u t )/At, where At is the time increment). Introducing Darcy's law with 
a permeability k assumed constant for simplicity, w = —kVp, the linearized 
formulation in matrix form is derived: 



K lp At K PP 



su \ r f-* \ _ r 

5p J \ AtF^* J \ AtF* nt 



(2.18) 



Each quantity (•) is linearised using the notations: (•) = (•) + S (•). Terms in 
(Eq. 2.18) are detailed here: 





/ {VN u } T DVN" + {VN"} T o-VN u ^ 

J V 


(2.19a) 




f {N P } T VN"& 

J V 


(2.19b) 




[ {VN p } T JcVN p * 

1 V 


(2.19c) 


T?int 
t u - 


1 {VN M } T <7 + {VN u } T N p p e A; 


(2.19d) 




f {VN P } T WN p p e + N p VN u v e dv 

J V 


(2.19e) 


■pext 


1 {N u } J tda 

J V 


(2.19f) 




j {N p } T WN p p e da 


(2.19g) 



where D is the matrix of material constants. In this work, the solid constituent is 
modelled with a Neo-Hooke model (see [4]): 

^NeoHooke = | (Ic ~ 3) - pL In J + ^ {In j f (2.20) 
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The formulation was implemented into an in-house finite element code. The 
observed aforementioned instabilities (presented in the Numerical example 
section) manifest in the form of spurious and non-physical oscillations in the 
pressure field in regions close to the free-flow boundaries. In the following 
section, we propose a stabilisation method to eliminate this issue. 



3 The GLS stabilization 

The spurious oscillations are stabilised using a Galerkin Least-Square (GLS) 
formulation, following [19]. This formulation was originally derived to solve 
geotechnical problems of fully and partially saturated soils. Although the 
derivation is extended here for finite deformations, the outcome is simpler owing 
to the fact that the permeability is assumed constant. A weighted least- square term 
R GLS , originating from the strong form of the fluid flow continuity equation (Eq. 
2.9) is derived and added onto the weak form (Eq. 2.17c): 

R = R u - AtR p + R GLS = (3.1) 

Starting off by defining the least-square term R GLS , where r* is a stabilisation 
factor that will subsequently be defined: 

R GLS = div (f + kVg^j T r* [div (v 5 + Wp)] dv (3.2) 

(Eq. 3.2) is rewritten as follows, taking into account the discretisation and the fact 
that the permeability is assumed constant and the pressure is linear: 

R GLS = | VN u f e } T r* {VN"u e } dv (3.3) 
Introducing the time integration scheme and defining t gls — t* / (At) 2 : 

R GLS = /{N"f e } T r GL5 {N M u e } dv (3.4) 

J V 

(Eq. 3.4) is linearised using the notations u e = u + 5u and writes in matrix form: 



K, 



+ K GLS 



AtK 



pp 



5u 
5p 



Fext 
u 

At F e p xt 



yint , y 



GLS-int 



At F* n * 



(3.5) 



8 



where: 



K^ s = f {VN M } T r GL5 {VN u } dv (3.6a) 

J V 

F GLS-int = f {VN"} T r GL5 {VN u u e } dv (3.6b) 

J v 

Finally, the stabilization factor is defined based on [19]: 

T GLS = ^ 

AkAt 

In [19], the parameter h is a length characterizing the element's size in the 
direction of the fluid flow. In the current application, it was simply taken as the 
radius of the element's circumsphere. As illustrated in the following section, this 
definition of r GLS stems from the fact that the solution needs greater stabilisation 
as the mesh gets coarser and as either or both the permeability and the size of the 
time-step decrease. In this work, due to the soft nature of the solid phase, it was 
not necessary to include a measure of material stiffness in the definition of t gls . 
Although this may be necessary for stiffer materials (e.g. [19] and [1]). 

It can be noted that the numerical integration is performed using a 4-point gaussian 
quadrature. The accuracy was verified with tests using up to 27 points. Reduced 
integration was also considered for the GLS terms, with no beneficial outcome. 



4 Numerical example 

The performance of the GLS stabilisation technique is assessed on a biphasic 
cylinder subjected to unconfined compression. With a 18mm radius, a thickness 
of 8mm and a solid phase defined with A = 0.2 MPa and /j, = 0.5 MPa, the cylinder 
can be thought of as an idealised human nucleus pulposus of the intervertebral 
disc when the permeability is set to k = 1 x 10 -3 mm 4 N -1 s -1 . 

The fluid flux q at the boundary is prescribed to zero on the vertical faces offering 
a lateral seal to the cylinder. The pressure p is set to zero on the top and bottom 
surfaces (see Fig. 2). These boundary conditions define the top and bottom 
surfaces as the only free-flow boundaries. This results in a near to uni-axial fluid 
flow through the depth of the cylinder. 

In what is herein referred to as the "reference loading", the top surface is displaced 
downwards at a rate of 2.5 fim.s^ 1 with time increments At = 6.4s, until the height 
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Figure 2: Loading and boundary conditions 



of the cylinder reduces by 1%. Such loading rates, together with permeabilities 
lower than 1 mm 4 N _1 s _1 , guaranties that the steady state is not reached instantly. 
In order to reduce the size of the problem, symmetry boundary conditions are 
applied onto a quarter cylinder. Analyses are undertaken on different meshes, the 
main characteristics of which are shown in Fig. 3. 

All results are plotted against nodal values gathered on the "reference line" defined 
in Fig. 2 and 3. 
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5261 5 


mesh 2 


16507 


11087 7 


mesh 3 


38424 


26622 9 


mesh 4 


78033 


55190 12 




Figure 3: Meshes characteristics (left) and mesh 2 (right) 



The low permeability hinders the fluid's ability to flow, defining two distinct 
regions associated with the load transfer mechanism. The first region, located 
near the free-flow boundaries, is dominated by solid deformation: the fluid does 
not have the ability to pressurize as it is squeezed out of the cylinder. This results 
in a lower fluid content in this region, explaining the peak strains observed near 
the top and bottom surfaces (see Fig. 4a). The second region, situated at the 
centre of the cylinder, is predominantly subjected to fluid pressurization (see Fig. 
4b) due to the fact that the low permeability is confining the fluid at the centre of 
the cylinder. 
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Cylinder Ordinate (mm) Cylinder Ordinate (mm) Cylinder Ordinate (mm) 

(a) Nodal strains (b) Nodal pressures (c) Influence of permeability 

Figure 4: Mesh 2: deformation mechanisms for k = 0.1 mm 4 N _1 s _1 (left and 
centre) and influence of permeability (right) 

As the permeability decreases, the boundary between solid- and 
pressure-dominated regions shifts towards the top and bottom surfaces (see 
Fig. 4c) and the level of pressurization rises. For a given mesh, when the 
permeability falls under a certain value (k < lxlO -1 mrn 4 N _1 s -1 for mesh 
2 in this example), the pressure profile starts to exhibit spurious oscillations 
near the free-flow boundaries (over 10% discrepancies for mesh 1 and mesh 2 
when k = 5x 10~ 2 mm 4 N _1 s _1 , and over 8% for all meshes when k = lx 10~ 2 
mm 4 N"V 1 ). Furthermore, the quality of the solution can be affected through 
the entire mesh as the near boundary oscillation propagates toward the centre for 
coarse meshes (e.g. mesh 1 in Fig. 5c). Finally, it was verified (in line with [20]) 
that decreasing the time-step exaggerates the pressure oscillations, although not 
presented here. 

Mesh refinement is the most natural and straightforward choice to overcome this 
issue, in particular when interested in accurately capturing the steep pressure 
gradients. As Fig. 5a and 5b illustrate, the non-physical pressure peaks can be 
removed by using meshes denser than that characterised by mesh 2. However, 
this is only a valid solution some cases. Fig. 5c illustrates that the spurious 
oscillations cannot always be reduced by reasonably- sized denser meshes for low 
permeability. 



11 




Figure 5: Effects of mesh refinement on standard Galerkin method 



The GLS stabilization offers substantial improvements to the solution. Fig. 6 
shows the benefits for four different meshes when k = 5x 10~ 3 mm 4 N _1 s _1 and k 
= 1 x 10~ 3 mm 4 N _1 s _1 . The primary enhancement is that all spurious oscillations 
observed in Fig. 5 have been stabilised, with the exception of mesh 2 where, 
when k = lxlO -3 mm 4 N' 1 s _1 , the discrepancies decreased from 27% to 6%. 
Additionally, when the accurate resolution of the near-boundary pressure gradient, 
defining the transition between the deformation- and pressure-driven regions, is 
not sought, the GLS formulation allows for coarser meshes to be used, since it 
also prevents the oscillations from propagating towards the centre. For example, 
when k = lxl0~ 3 mm 4 N~ 1 s~ 1 , mesh 1 (8000 nodes) with GLS offers similar 
performances to mesh 3 (38000) without the stabilisation (compare Fig. 5c and 
6b). Finally, it is important to notice that the GLS stabilisation is only having 
a damping effect on the spurious oscillations, while leaving stable solutions 
unaffected (compare Fig. 5b and 6a). 




Figure 6: Effects of mesh refinement on GLS stabilisation 



A few observations can be made to support the choice of the stabilization factor 
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t in (Eq. 3.7). First, Fig. 6 highlights the fact that the stabilization 
performs equally well irrespective of the change of mesh and permeability, 
giving confidence in the way the element's characteristic size and permeability 
are accounted for. The impact of the time- step on the stabilization was also 
investigated: it was verified that if r GLS is not inversely proportional to At, 
stabilisation is not possible. Also, in simulations not presented here, it was 
confirmed that the size of the time-step (At = {0.64s, 3.2s, 6.4s, 8s, 32s}) does 
not affect the quality of the stabilised solution. Finally, it was also verified 
(again not shown here) that changing the loading rate (1.25//m.s _1 , 2.5 fim.s" 1 , 
6.25fim.s~ 1 ) did not affect the degree of peak pressure oscillations for the 
stabilised results. 

Performance of the GLS stabilization was initially assessed for greater levels of 
deformation. As Fig. 7 shows, the level of pressure discrepancies reduces as the 
compressive strain increases, which is in line with the findings in [19] and [1], 
where oscillations are reported to occur at the "early stage" of the consolidation 
problem. Although the GLS stabilisation also performs well at higher strains 
(see Fig. 7b), this observation motivated the choice to present results at 1% 
compression throughout this numerical example. 




-4 -3 -2 -1 1 2 3 4 -4 -3 -2 -1 1 2 3 4 



Cylinder Ordinate (mm) Cylinder Ordinate (mm) 

(a) Standard Galerkin (b) GLS 

Figure 7: Effects of higher strain on stability for mesh 2 

A sensitivity study was performed to charaterise the parameter h in (Eq. 
3.7). Several combinations of the radius of the circumsphere, the shortest and 
longest edges of the tetrahedron were tested without noticeable and consistent 
improvement to the overall solution. 
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5 Conclusion 



It was observed that an hyperelastic biphasic model, implemented in a finite 
element framework with Taylor-Hood tetrahedral elements, exhibits non-physical 
pressure oscillations for low permeabilities. A Galerkin least-square formulation 
was derived for finite deformations in order to stabilise these oscillations. 

In the context of constant permeability and near to uni-axial fluid flow, the current 
formulation shows good results. It eliminates the spurious oscillations for most 
meshes (and damp the oscillations for others meshes) and also prevents these 
oscillations from propagating towards the centre of the medium as reported for 
very coarse meshes. The solution scheme proved to be robust when tested against 
various mesh densities, permeabilities, loading rates, compressive strains and time 
steps. It is also worth mentioning that the benefits of this formulation come at 
minimal computational cost, as no additional degrees of freedom are required. 

For more complex fluid flow situations, information regarding the fluid 
flow directionality should be included in the derivation of the element's 
characteristic size. In [19], this has been considered for 2D and for 3D brick 
elements. It will prove more challenging for 3D tetrahedral elements. For 
extension to strain-dependent permeability, derivation of the stabilisation terms 
is straightforward but would result in a non-symmetrical system that would be 
computationally more expensive to solve. 
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